function [rr_sqn,varargout] = Path_ctrf(mu_init,Q_init,varargin)

Par = varargin{1};
Step  = varargin{2};
Trans = varargin{3};

steps      = Step{1};
Trans_innov = Step{2};
Trans_fail = Step{3};
Trans_exo  = Step{4};

Param = structvars(Par);
for ctr0 = 1:size(Param,1); eval(Param(ctr0,:)); end

Transit = structvars(Trans);
for ctr0 = 1:size(Transit,1); eval(Transit(ctr0,:)); end

tic
%============== Steady State =====================

QA = [Q_init zeros(2*mbar+1,time_ss-1)];
QB = diag(lmd.^steps)*flip(QA,1);

% ---- Solve for SS Values ----

rr_init = log([.02730; .02730; 1; 1]);

if exist('rr_guess.mat','file')>0

    try load('rr_guess.mat'); rr_init = log(rr_ss); catch; end

end    
      
opt_steady = optimset('Display','final');
        
steady_new = @(N)Steady(N,Par,Step,Trans,QA(:,1));
[rr_ss,fval_ss,exit_ss,output_ss] = fsolve(steady_new,rr_init,opt_steady);

% keyboard

if exit_ss~=1 && max(abs(fval_ss))>2.5*10^-3
    
    disp(' ')
    disp(['          A 2nd Round of R-loop !!! Starting at ' ...
            datestr(now, 'HH:MM:SS') ' !!!']);
    disp(' ')
   
    rr_init = real(rr_ss);
    opt_steady = optimset('Display','final','TolX',10^-10);
    
    steady_new = @(N)Steady(N,Par,Step,Trans,QA(:,1));
    [rr_ss,fval_ss,exit_ss,output_ss] = fsolve(steady_new,rr_init,opt_steady);
    
end

if exit_ss~=1 && max(abs(fval_ss))>2.5*10^-3
    
    disp(' ')
    disp(['          A 3rd Round of R-loop !!! Starting at ' ...
            datestr(now, 'HH:MM:SS') ' !!!']);
    disp(' ')
   
    rr_init = log([.008; .008; 1; 1]);
    opt_steady = optimset('Display','final','TolX',10^-10);
    
    steady_new = @(N)Steady(N,Par,Step,Trans,QA(:,1));
    [rr_ss,fval_ss,exit_ss,output_ss] = fsolve(steady_new,rr_init,opt_steady);
    
end

if ~isreal(rr_ss)
    
    rr_sqn = 0;
    varargout = cell(max(nargout,1)-1,1);
    
    disp(' ')
    disp('          R in SS Problem - Returning !!!');
    disp(' ')
    
    return
    
end

[~, set_ss, prwg_ss] = Steady(rr_ss,Par,Step,Trans,QA(:,1));

rr_ss   = exp(rr_ss);

if exit_ss==1; save rr_guess rr_ss; end

if exist('rr_diff.mat','file')>0
    
    delete('rr_diff.mat')
    
end

disp(' ')
disp(['          R in SS solved !!! Finishing at ' ...
        datestr(now, 'HH:MM:SS') ' !!!']);
disp(' ')
disp('         FINAL      ')
disp('     rr_ss    fval_ss')
disp(' ')
disp([rr_ss fval_ss])
disp(' ')
toc_ss = toc;

%============== Transition Path =====================

if exist('RQ.mat','file')>0
    
    load RQ.mat;

    rr_sqn = [RQ(1:2,1:end-1) rr_ss(1:2)];
    qq_sqn = [RQ(3:4,1:end-1) rr_ss(3:4)];
    
else

    rr_sqn = rr_ss(1:2)*ones(1,time_ss);
    qq_sqn = rr_ss(3:4)*ones(1,time_ss);
    
end

% rr_sqn = rr_ss(1:2)*ones(1,time_ss);
% qq_sqn = rr_ss(3:4)*ones(1,time_ss);

V_iA = set_ss(:,1)*[zeros(1,time_ss-1) 1];
V_iB = set_ss(:,2)*[zeros(1,time_ss-1) 1];

X_iA = set_ss(:,3)*[zeros(1,time_ss-1) 1];
X_iB = set_ss(:,4)*[zeros(1,time_ss-1) 1];
X_eA = set_ss(:,5)*[zeros(1,time_ss-1) 1];
X_eB = set_ss(:,6)*[zeros(1,time_ss-1) 1];

mu_share = zeros(length(mu_init),time_ss);
mu_share(:,1) = mu_init;

prft_fnc_A = prwg_ss(:,1)*[zeros(1,time_ss-1) 1];
prft_fnc_B = prwg_ss(:,2)*[zeros(1,time_ss-1) 1];
wage_dom_A = prwg_ss(:,3)*[zeros(1,time_ss-1) 1];
wage_dom_B = prwg_ss(:,4)*[zeros(1,time_ss-1) 1];
                
wq_sqn = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(1,:).^(-Par.bet);
            (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn(2,:).^(-Par.bet)];  

% ===== Cutoff series =====

Qrat  = qq_sqn(1,:)./qq_sqn(2,:);

mrealX = (1-bet)*(+1/bet*log((1+Par.kapa)*(1+Par.trf_B))-log(Qrat))/log(lmd); 
mrealM = (1-bet)*(-1/bet*log((1+Par.kapa)*(1+Par.trf_A))-log(Qrat))/log(lmd); 

madjX = round(mrealX)-1; % largest step with zero exports in A
madjX = min(max(madjX,-(mbar+1)),mbar);

madjM = -(round(mrealM)+1); % smallest step with zero imports in A (absolute val)
madjM = min(max(madjM,-(mbar+1)),mbar);

mcutX = madjX+1;
mcutM = madjM+1;

epsX = 1-(mrealX-round(mrealX)+0.5);
epsM = -(mrealM-round(mrealM)-0.5); 

Tol_rr = 10^-6;
iter_r = 1;
iter_N = 800;
diff_q = 1;
diff = 1;
load Tol_rr

while (diff > Tol_rr || diff_q > Tol_rr) && iter_r < iter_N
    iter_r

    vv_iA_next = set_ss(:,1);
    vv_iB_next = set_ss(:,2);

    % ===== Retrieve R&D decisions X =====

    for ctrx = time_ss-1:-1:1

        % ===== Profit functions =====

        prft_A = [ones(Par.mbar-mcutX(ctrx),1)*(Par.pi_dom*L_A+Par.pi_frn*L_B);
                ones((Par.mbar-madjX(ctrx))>0,1)'*Par.pi_dom*Par.L_A+epsX(ctrx)*Par.pi_frn*Par.L_B;
                ones(madjX(ctrx)+madjM(ctrx)+1,1)*Par.pi_dom*L_A; 
                ones((Par.mbar-madjM(ctrx))>0,1)'*epsM(ctrx)*Par.pi_dom*Par.L_A;
                zeros(Par.mbar-mcutM(ctrx),1)];

        prft_B = [ones(Par.mbar-mcutM(ctrx),1)*(Par.pi_dom*L_B+Par.pi_frn_toUS*L_A);
                ones((Par.mbar-madjM(ctrx))>0,1)'*Par.pi_dom*Par.L_B+(1-epsM(ctrx))*Par.pi_frn_toUS*Par.L_A;
                ones(madjX(ctrx)+madjM(ctrx)+1,1)*Par.pi_dom*L_B; 
                ones((Par.mbar-madjX(ctrx))>0,1)'*(1-epsX(ctrx))*Par.pi_dom*Par.L_B;
                zeros(Par.mbar-mcutX(ctrx),1)];

        prft_fnc_A(:,ctrx) = prft_A;
        prft_fnc_B(:,ctrx) = prft_B;        

        wage_dom_A(:,ctrx) = [ones(1,Par.mbar+1+madjM(ctrx)) ones((Par.mbar-madjM(ctrx))>0,1)*epsM(ctrx) ...
                                zeros(1,Par.mbar-mcutM(ctrx))]'; 

        wage_dom_B(:,ctrx) = [ones(1,Par.mbar+1+madjX(ctrx)) ones((Par.mbar-madjX(ctrx))>0,1)*(1-epsX(ctrx)) ...
                                zeros(1,Par.mbar-mcutX(ctrx))]'; 

        % ===== R&D Decisions using Vt+dt =====

        X_iA(:,ctrx) = ((Trans_innov'*vv_iA_next-...
                        vv_iA_next)/(alfa_A*(1-tau_A))).^(1/(gama_A-1));
        X_iB(:,ctrx) = ((Trans_innov'*vv_iB_next-...
                        vv_iB_next)/(alfa_B*(1-tau_B))).^(1/(gama_B-1));

        X_eA(:,ctrx) = ((Trans_innov'*vv_iA_next...
                        -0)/alfa_eA).^(1/(gama_eA-1));
        X_eB(:,ctrx) = ((Trans_innov'*vv_iB_next...
                        -0)/alfa_eB).^(1/(gama_eB-1));

        % ===== Value Functions at t =====


        v_iA = (prft_fnc_A(:,ctrx)*wq_sqn(1,ctrx)^(1-1/bet)-...
               (1-tau_A)*alfa_A.*X_iA(:,ctrx).^gama_A/gama_A+...
               X_iA(:,ctrx).*(Trans_innov'*vv_iA_next-...
               vv_iA_next)-X_eA(:,ctrx).*vv_iA_next+...
               flip(X_eB(:,ctrx)+X_iB(:,ctrx)).*...
               (Trans_fail*vv_iA_next-vv_iA_next)+...
               dlt*Trans_exo*vv_iA_next-...
               (rr_sqn(1,ctrx)-1/dt)*vv_iA_next)*dt;

        v_iB = (prft_fnc_B(:,ctrx)*wq_sqn(2,ctrx)^(1-1/bet)-...
               (1-tau_B)*alfa_B.*X_iB(:,ctrx).^gama_B/gama_B+...
               X_iB(:,ctrx).*(Trans_innov'*vv_iB_next-...
               vv_iB_next)-X_eB(:,ctrx).*vv_iB_next+...
               flip(X_eA(:,ctrx)+X_iA(:,ctrx)).*...
               (Trans_fail*vv_iB_next-vv_iB_next)+...
               dlt*Trans_exo*vv_iB_next-...
               (rr_sqn(2,ctrx)-1/dt)*vv_iB_next)*dt;

        % ===== Update Vt+dt moving to previous period =====

        vv_iA_next = v_iA;
        vv_iB_next = v_iB;

        V_iA(:,ctrx) = v_iA;
        V_iB(:,ctrx) = v_iB;

    end

    % ===== Generate Path Using X =====    

    for ctrp = 1:time_ss-1

        % ===== Transition Matrices =====

        Qtr_A_idom = Qloc_idom*diag(X_iA(:,ctrp));
        Qtr_A_ifrn = Qloc_ifrn_L*(diag(X_iB(:,ctrp))*Qloc_frn_R);
        Qtr_A_edom = Qloc_edom*diag(X_eA(:,ctrp));
        Qtr_A_efrn = Qloc_efrn_L*(diag(X_eB(:,ctrp))*Qloc_frn_R);

        Qtr_B_idom = Qloc_idom*diag(X_iB(:,ctrp));
        Qtr_B_ifrn = Qloc_ifrn_L*(diag(X_iA(:,ctrp))*Qloc_frn_R);
        Qtr_B_edom = Qloc_edom*diag(X_eB(:,ctrp));
        Qtr_B_efrn = Qloc_efrn_L*(diag(X_eA(:,ctrp))*Qloc_frn_R);

        Trans_matA = (Qtr_A_idom+Qtr_A_ifrn+Qtr_A_edom+Qtr_A_efrn+Qtr_exo)*dt;
        Trans_matB = (Qtr_B_idom+Qtr_B_ifrn+Qtr_B_edom+Qtr_B_efrn+Qtr_exo)*dt;

        QA(:,ctrp+1) = QA(:,ctrp)+Trans_matA*QA(:,ctrp);
        QB(:,ctrp+1) = QB(:,ctrp)+Trans_matB*QB(:,ctrp);

    end 

    QwA = sum(wage_dom_A.*QA+flip(1-wage_dom_B).*QA*(1+kapa)/((1+kapa)*(1+trf_B))^(1/bet));
    QwB = sum(wage_dom_B.*QB+flip(1-wage_dom_A).*QB*(1+kapa)/((1+kapa)*(1+trf_A))^(1/bet));

    qq_sqn1 = [sum(QA)./QwA; sum(QB)./QwB];
                
    wq_sqn1 = [(1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn1(1,:).^(-Par.bet);
                 (1-Par.bet)*Par.eta^(Par.bet-1)*qq_sqn1(2,:).^(-Par.bet)];  

    xk_sqn1 = [(Par.pi_frn/Par.bet)^(1/(1-Par.bet))./wq_sqn1(1,:).^(1/Par.bet)*Par.L_B;
             (Par.pi_frn_toUS/Par.bet)^(1/(1-Par.bet))./wq_sqn1(2,:).^(1/Par.bet)*Par.L_A];

    IA = sum(prft_fnc_A.*QA).*wq_sqn1(1,:).^(1-1/bet)-...
        sum(alfa_A/gama_A*X_iA.^gama_A.*QA)-...
        sum(alfa_eA/gama_eA*X_eA.^gama_eA.*QA)+...
        pi_dom/(1-bet)*L_A*wq_sqn1(1,:).^(1-1/bet).*sum(wage_dom_A.*QA)+...
        pi_frn_toUS/(1-bet)*L_A*wq_sqn1(2,:).^(1-1/bet).*sum(flip(1-wage_dom_A).*QB)+...
        wq_sqn1(1,:).*sum(QA)+...
        Par.rev_US*trf_A*((1+trf_A)*(1+kapa)*eta/(1-bet))*wq_sqn1(2,:).*xk_sqn1(2,:).*sum(flip(1-wage_dom_A).*QB);  
    IB = sum(prft_fnc_B.*QB).*wq_sqn1(2,:).^(1-1/bet)-...
        sum(alfa_B/gama_B*X_iB.^gama_B.*QB)-...
        sum(alfa_eB/gama_eB*X_eB.^gama_eB.*QB)+...
        pi_dom/(1-bet)*L_B*wq_sqn1(2,:).^(1-1/bet).*sum(wage_dom_B.*QB)+...
        pi_frn/(1-bet)*L_B*wq_sqn1(1,:).^(1-1/bet).*sum(flip(1-wage_dom_B).*QA)+...
        wq_sqn1(2,:).*sum(QB)+...
        Par.rev_FN*trf_B*((1+trf_B)*(1+kapa)*eta/(1-bet))*wq_sqn1(1,:).*xk_sqn1(1,:).*sum(flip(1-wage_dom_B).*QA); 

    g_t = [(IA(2:end)-IA(1:end-1))./IA(1:end-1); ...
           (IB(2:end)-IB(1:end-1))./IB(1:end-1)];
       
    % Finishing steps

    diff = max(max(abs(psy/dt*g_t(:,1:end)+ro-rr_sqn(:,1:end-1))))
    diff_q = max(max(abs([sum(QA)./QwA; sum(QB)./QwB]-qq_sqn)))

    rr_sqn1 = rr_sqn;
    rr_sqn  = [num_r*(psy/dt*g_t+ro)+(1-num_r)*rr_sqn(:,1:end-1) rr_sqn(:,end)];
    qq_sqn  = qq_sqn+num_r*(qq_sqn1-qq_sqn);

    Qrat  = qq_sqn(1,:)./qq_sqn(2,:);

    madjX0 = madjX;
    madjM0 = madjM;

    mrealX = (1-bet)*(+1/bet*log((1+kapa)*(1+trf_B))-log(Qrat))/log(lmd); 
    mrealM = (1-bet)*(-1/bet*log((1+kapa)*(1+trf_A))-log(Qrat))/log(lmd); 

    madjX = round(mrealX)-1; % largest step with zero exports in A
    madjX = min(max(madjX,-(mbar+1)),mbar);

    madjM = -(round(mrealM)+1); % smallest step with zero imports in A (absolute val)
    madjM = min(max(madjM,-(mbar+1)),mbar);

    mcutX = madjX+1;
    mcutM = madjM+1;

    epsX = 1-(mrealX-round(mrealX)+0.5);
    epsM = -(mrealM-round(mrealM)-0.5);  

    diff_xm = (max(abs(madjX-madjX0)+abs(madjM-madjM0))>0)

    iter_r = iter_r+1;

    if diff_q>10; iter_r=iter_N; end

end
        
% ===== Leadership Shares =====

Tmu_A = 0;

for ctrm = burn_ss+1:time_ss
        
    % ===== Transition Matrices =====
    
    Mutr_A_idom = Muloc_idom*diag(X_iA(:,ctrm));
    Mutr_A_ifrn = Muloc_ifrn_L*(diag(X_iB(:,ctrm))*Muloc_frn_R);
    Mutr_A_edom = Muloc_edom*diag(X_eA(:,ctrm));
    Mutr_A_efrn = Muloc_efrn_L*(diag(X_eB(:,ctrm))*Muloc_frn_R);
    
    Mutr_B_idom = Muloc_idom*diag(X_iB(:,ctrm));
    Mutr_B_ifrn = Muloc_ifrn_L*(diag(X_iA(:,ctrm))*Muloc_frn_R);
    Mutr_B_edom = Muloc_edom*diag(X_eB(:,ctrm));
    Mutr_B_efrn = Muloc_efrn_L*(diag(X_eA(:,ctrm))*Muloc_frn_R);
    
    Trans_muA = (Mutr_A_idom+Mutr_A_ifrn+Mutr_A_edom+Mutr_A_efrn+Mutr_exo)*dt;
    Trans_muB = (Mutr_B_idom+Mutr_B_ifrn+Mutr_B_edom+Mutr_B_efrn+Mutr_exo)*dt;
    
    Tmu_A = Tmu_A+Trans_muA;
    
    mu_share(:,ctrm+1) = mu_share(:,ctrm)+Trans_muA*mu_share(:,ctrm);

end
toc_tr = toc;

%============== Finishing =====================

ind_neg = 0;
if sum(sum((X_iA<0)+(X_iB<0)+(X_eA<0)+(X_eB<0)))>0
    
    ind_neg = 1;
    
end

if sum(sum(rr_sqn<0))>0 || sum(sum(isnan(rr_sqn)))>0 || sum(sum(1-isreal(rr_sqn)))>0

    ind_neg = 2;

end

if diff_xm

    ind_neg = 3;

end

if iter_r == iter_N

    ind_neg = 4;

end

if diff_q > 10

    ind_neg = 5;

end

%=============== Add Span to the solution already at SS =================
% Fix r, Xe at SS, continue from last period Q's found above for "time_add" long! 

rr_sqn = [rr_sqn rr_ss(1:2)*ones(1,time_add)];
l_rr   = length(rr_sqn);

QA_long = [QA zeros(size(QA,1),l_rr-time_ss)];
QB_long = [QB zeros(size(QA,1),l_rr-time_ss)];  
mu_long = [mu_share zeros(size(mu_share,1),l_rr-(time_ss+1))];

if if_long

    Qtr_A_idom = Qloc_idom*diag(X_iA(:,end));
    Qtr_A_ifrn = Qloc_ifrn_L*(diag(X_iB(:,end))*Qloc_frn_R);
    Qtr_A_edom = Qloc_edom*diag(X_eA(:,end));
    Qtr_A_efrn = Qloc_efrn_L*(diag(X_eB(:,end))*Qloc_frn_R);

    Qtr_B_idom = Qloc_idom*diag(X_iB(:,end));
    Qtr_B_ifrn = Qloc_ifrn_L*(diag(X_iA(:,end))*Qloc_frn_R);
    Qtr_B_edom = Qloc_edom*diag(X_eB(:,end));
    Qtr_B_efrn = Qloc_efrn_L*(diag(X_eA(:,end))*Qloc_frn_R);

    Trans_matA = (Qtr_A_idom+Qtr_A_ifrn+Qtr_A_edom+Qtr_A_efrn+Qtr_exo)*dt;
    Trans_matB = (Qtr_B_idom+Qtr_B_ifrn+Qtr_B_edom+Qtr_B_efrn+Qtr_exo)*dt;
    %

    Mutr_A_idom = Muloc_idom*diag(X_iA(:,end));
    Mutr_A_ifrn = Muloc_ifrn_L*(diag(X_iB(:,end))*Muloc_frn_R);
    Mutr_A_edom = Muloc_edom*diag(X_eA(:,end));
    Mutr_A_efrn = Muloc_efrn_L*(diag(X_eB(:,end))*Muloc_frn_R);

    Trans_muA = (Mutr_A_idom+Mutr_A_ifrn+Mutr_A_edom+Mutr_A_efrn+Mutr_exo)*dt;

    for ctradd = time_ss:l_rr-1   

        QA_long(:,ctradd+1) = QA_long(:,ctradd)+Trans_matA*QA_long(:,ctradd);
        QB_long(:,ctradd+1) = QB_long(:,ctradd)+Trans_matB*QB_long(:,ctradd);

        mu_long(:,ctradd+1) = mu_long(:,ctradd)+Trans_muA*mu_long(:,ctradd);

    end

end

RQ = [rr_sqn; qq_sqn];
% save RQ RQ

varargout{1} = [X_iA; X_iB];
varargout{2} = [V_iA; V_iB];
varargout{3} = [X_eA; X_eB];
varargout{4} = [QA; QB];
varargout{5} = [QA_long; QB_long];
varargout{6} = mu_long;
varargout{7} = ind_neg;
varargout{8} = [prft_fnc_A; prft_fnc_B];
varargout{9} = [wage_dom_A; wage_dom_B];
varargout{10} = [mcutX; mcutM];